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Abstract 

The concept of a universal algorithm is discussed. Examples of this 
kind of algorithms are presented. Software implementations of such 
algorithms in C ++ -type languages are discussed together with means 
that provide for computations with an arbitrary accuracy. Particular 
emphasis is placed on universal algorithms of linear algebra over semir- 
ings. 



INTRODUCTION 

Modern achievements in software development and mathematics make us 
consider numerical algorithms and their classification from a new point of 
view. Conventional numerical algorithms are oriented to software (or hard- 
ware) implementation with the use of the floating point arithmetic and fixed 
accuracy. However, it is often desirable to perform computations with vari- 
able (and arbitrary) accuracy. For this purpose, algorithms are required 
that are independent of the accuracy of computations and of a particular 
computer representation of numbers. In fact, many algorithms are not only 
independent of the computer representation of numbers, but also of concrete 
mathematical (algebraic) operations on data. In this case, operations may 
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be considered as variables. Such algorithms are implemented by generic pro- 
grams based on the abstract data types technique (abstract data types are 
defined by the user, in addition to predefined types of the language used). 
The corresponding program tools appeared as early as in Simula-67, but mod- 
ern object-oriented languages (like C ++ , see, e.g. |], @]) are more convenient 
for generic programming. 

The concept of a generic program was introduced by many authors; for 
example, in ||, such programs were called program schemes. In this pa- 
per, we discuss universal algorithms implemented as generic programs and 
their specific features. This paper is closely related to papers f§, K, in which 
the concept of a universal algorithm was defined and software and hardware 
implementation of such algorithms was discussed in connection with prob- 
lems of idempotent mathematics (|, |]. In this paper, the emphasis is placed 
on software implementation of universal algorithms, computations with ar- 
bitrary accuracy, universal algorithms of linear algebra over semirings, and 
their implementation in C ++ . 

1. UNIVERSAL ALGORITHMS 

Computational algorithms are constructed on the basis of certain ba- 
sic operations. Basic operations manipulate data that describe "numbers". 
These "numbers" are elements of a "numerical domain", i.e., a mathemati- 
cal object like the field of real numbers, the ring of integers, or an idempo- 
tent semiring of numbers (idempotent semirings and their role in idempotent 
mathematics are discussed in |4], || and below in this paper) . In every partic- 
ular computation, elements of the numerical domains are replaced by their 
computer representations, i.e., by elements of certain finite models of these 
domains. Examples of models that can be conveniently used for computer 
representation of real numbers are provided by various modifications of float- 
ing point arithmetics, approximate arithmetics of rational numbers |7|], and 
interval arithmetics. The difference between mathematical objects ("ideal" 
numbers) and their finite models (computer representations) results in com- 
putational (e.g., rounding) errors. 
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An algorithm is called universal if it is independent of a particular nu- 
merical domain and (or) of its computer representation. A typical exam- 
ple of a universal algorithm is computation of the scalar product (x, y) of 
two vectors x = (xx, . . . , x n ) and y = (y\, . . . , y n ) by the formula (x, y) = 
x\y\ + . . . + x n y n . This algorithm (formula) is independent of a particular 
domain and its computer implementation, since the formula is defined for 
any semiring. It is clear that one algorithm can be more universal than 
another. For example, the simplest rectangular formula provides the most 
universal algorithm for numerical integration; indeed, this formula is valid 
even for idempotent integration (over any idempotent semiring 0). Other 
quadrature formulas (e.g., combined trapezoid or Simpson formulas) are inde- 
pendent of computer arithmetics and can be used (e.g., in the iterative form) 
for computations with arbitrary accuracy. In contrast, algorithms based on 
Gauss- Jacobi formulas are designed for fixed accuracy computations: they 
include constants (coefficients and nodes of these formulas) defined with fixed 
accuracy. Certainly, algorithms of this type can be made more universal by 
including procedures for computing the constants; however, this results in an 
unjustified complication of the algorithms. 

Computer algebra algorithms used in such systems as Mathematica, Maple, 
REDUCE, and others are higly universal. Standard algorithms used in lin- 
ear algebra can be rewritten in such a way that they will be valid over any 
field and complete idempotent semiring (including semirings of intervals; see 
[§, |§, where an interval version of the idempotent linear algebra and the 
corresponding universal algorithms are discussed). 

As a rule, iterative algorithms (beginning with the successive approxi- 
mation method) for solving differential equations (e.g., methods of Euler, 
Euler-Cauchy, Runge-Kutta, Adams, a number of important versions of the 
difference approximation method, and the like), methods for calculating ele- 
mentary and some special functions based on the expansion in Taylor's series 
and continuous fractions (Pade approximations) and others are independent 
of the computer representation of numbers. 
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2. UNIVERSAL ALGORITHMS AND ACCURACY 
OF COMPUTATIONS 



Calculations on computers usually are based on a floating-point arith- 
metic with a mantissa of a fixed length; i.e., computations are performed 
with fixed accuracy. Broadly speaking, with this approach only the rela- 
tive rounding error is fixed, which can lead to a drastic loss of accuracy 
and invalid results (e.g., when summing series and subtracting close num- 
bers). On the other hand, this approach provides rather a high speed of 
computations. Many important numerical algorithms are designed to use 
a floating-point arithmetic (with fixed accuracy) and ensure the maximum 
computation speed. However, these algorithms are not universal. The above 
mentioned Gauss- Jacobi quadrature formulas, computation of elementary 
and special functions on the basis of the best polynomial or rational approxi- 
mations or Pade-Chebyshev approximations, and some others belong to this 
type. Such algorithms use nontrivial constants specified with fixed accuracy. 

Recently, problems of accuracy, reliability, and authenticity of computa- 
tions (including the effect of rounding errors) have come to the fore; in part, 
this fact is related to the ever-increasing performance of computer hardware. 
When errors in initial data and rounding errors strongly affect the computa- 
tion results (ill-posed problems, analysis of stability of solutions, etc.), it is 
often useful to perform computations with improved and variable accuracy. 
In particular, the rational arithmetic, in which the rounding error is speci- 
fied by the user , can be used for this purpose. This arithmetic is a useful 
complement to the interval analysis flO ]. The corresponding computational 
algorithms must be universal (in the sense that they must be independent of 
the computer representation of numbers). 



4. MATHEMATICS OF SEMIRINGS 



A broad class of universal algorithms is related to the concept of a semir- 
ing. We reiterate here the definition of a semiring (see, e.g., |TI[). Let 
S be a set on which associative binary operations © and 0, called addi- 
tion and multiplication, respectively, are defined. We assume that addition 
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is commutative and that multiplication is distributive over addition; i.e., 
x (y © z) = (x y) © (x © z) and (x © y) z = (x z) © (y *) for all 
x,y,z G S 1 . In this case, 5 is called a semiring. We assume that the semiring 
S contains identity 1 and zero 0; i.e., 10i = i01 = i and © x = x, 
001 = 100 = 0; in addition, 0^1. As is customary, we sometimes omit 
the multiplication symbol. 

A semiring S is called commutative if multiplication © is commutative. 
A semiring S is called idempotent if x © x = x for all x 6 S. If a semiring is 
a group under addition, it is called a ring (in this case, it cannot be idempo- 
tent). If every nonzero element of a commutative ring (semiring) is invertible 
under multiplication, this ring (semiring) is called a field (semifield) . 

The best known and most important examples of semirings are "numer- 
ical" semirings consisting of real numbers. For example, the set R of all 
real numbers is a field under ordinary arithmetic operations; i.e., © = +, 
= ■, = 0, 1 = 1. The set R max = R U {— oo} equipped with operations 
© = max and = + provides an example of an idempotent semiring (and 
semifield). Here = — oo and 1 = 0. This semifield is often called the Max- 
Plus algebra. The semiring R m i n = R U {+00} equipped with operations 
© = min and = + is isomorphic to the Max-Plus algebra. Here = +00 
and 1 = 0. Another example is the set stiax,min consisting of the elements 
of an interval [a, b], where —00 < a < b < +00, equipped with operations 
© = max and © = min; here = a and 1 = b. This commutative semiring 
is not a semifield. 

An important example of a noncommutative semiring is the set Mat n (<S') 
of all matrices of order n x n with elements from a commutative semiring 
S with ordinary standard operations. The sum of matrices A = (a^) and 
B = {pij) is the matrix A © B = (a^- © bij), and the product of these matrices 
is the matrix AB = (©£ =1 © 6^), where i,j = l,...,n. Operations 
on rectangular matrices can be defined similarly. Zero O and identity / in 
M&t n (S) are defined in the conventional way. If the semiring S is idempotent, 
then Mat n (5') is also idempotent. Many other important examples can be 
found in § - 1, [|, 1, [0. 

On any idempotent semiring, a canonical partial order 4 is defined by the 
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following rule: x =4. y is equivalent to x © y = y. Moreover, x © y = sup{x, y} 
with respect to the canonical order. The canonical order is compatible with 
the semiring addition and multiplication in the common way. For the semir- 
ings R max and S^^l min , the canonical order coincides with the standard order 
< defined on the set of real numbers; for the semiring R m i n , it is inverse to 
the standard order. 

There exists a (heuristic) correspondence between important, useful, and 
interesting constructs and results of traditional mathematics over fields and 
similar constructs and results of idempotent mathematics (i.e., mathemat- 
ics over idempotent semirings). This idempotent correspondence principle 
is closely related to the Bohr correspondence principle in quantum mechan- 
ics. Traditional mathematics can be considered as a "quantum" theory and 
idempotent mathematics as its "classical" analogue (see ||). Consistent ap- 
plication of the idempotent correspondence principle leads to various and 
surprising results, including a methodology for constructing universal algo- 
rithms and patenting computer devices ||, ||. 

The fundamental equations in quantum theory are linear (superposition 
principle). There is an idempotent version of the superposition principle ||): 
the Hamilton- Jacobi equation, i.e., the basic (nonlinear) equation of classical 
mechanics, can be considered as linear over the semiring R m i n ; various mod- 
ifications of the Bellman equation, i.e., the basic equation of optimization 
theory, are also linear over appropriate idempotent semirings. For example, 
the finite-dimensional time-independent Bellman equation can be written as 

X = AQX®B, (1) 

where A is a square matrix with elements from an idempotent semiring 5* 
and X and B are vectors (or matrices) with elements from S. The solution 
X is found from (1) when A and B are given. 

In particular, standard problems in dynamic programming correspond to 
the case R max , and the well-known shortest path problem corresponds to 
S = Rmax- It is shown in [|l^] that the principal optimization algorithms for 
finite graphs correspond to standard methods for solving systems of linear 
equations of form (1) over semirings. The Bellman algorithm for the shortest 
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path problem corresponds to a semiring version of the Jacobi method; the 
Ford algorithm corresponds to the Gauss-Seidel iterative method; and so 
on. These algorithms are universal and may be used for solving linear al- 
gebra problems over a broad class of semirings that includes all idempotent 
semirings and all fields. 

Idempotent analogues of standard numerical algorithms are very impor- 
tant and can be used systematically for solving, for example, optimization 
problems. Linear algebra algorithms are of prime importance, since standard 
infinite-dimensional linear problems over semirings can be reduced to finite- 
dimensional (or finite) approximations, and nonlinear algorithms can often 
be approximated by linear ones. 

We note that the available methods used for parallelizing linear algebra 
algorithms can be applied to their semiring analogues. 



The most important linear algebra problem is solving the system of linear 
equations 



where A is a matrix with elements from the basic field and X and B are 
vectors (or matrices) with elements from the same field. It is required to find 
X if A and B are given. If A in (2) is not the identity matrix J, then system 
(2) can be written in form (1), i.e., 



It is well known that form (1) or (1') is convenient for using the successive 
approximation method. Applying this method with the initial approximation 
X = 0, we obtain the solution 



5. UNIVERSAL LINEAR ALGEBRA 
ALGORITHMS OVER SEMIRINGS 



AX = B, 



(2) 



X = AX + B. 



(!') 



X = A*B, 



(3) 



where 



A* = I + A + A 2 



+ ... + A 



+ ... 



(4) 
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On the other hand, it is clear that 



A* = (I-A)-\ (5) 

if the matrix I — A is invertible. The inverse matrix (/ — A)^ 1 can be 
considered clS 8b TQ gularized sum of the formal series (4). 

The above considerations can be extended to a broad class of semirings. 
The unary operation A i— > A* in Mat n (S') is defined (partially) if a unary 
(partial) operation x t— > x*, called closure, is defined on the semiring S such 
that the identity 

x* = l©(x*Ox) = l©(x©x*) (6) 
holds true if x* is defined. It follows from (6) that 

™* — 1 £d t (X\ en en <r n fft ^*^. n +i 

X — ± W X KiJ X KD . . . KD X Kiy X X 

for any positive integer n; thus, x* can be considered as a regularized sum of 
the formal series 

™* 1 (T\ rf, (T\ rJ^ (T\ (T\ rf n (T\ 

X — ± KD X KD X KD . . . KD X KD . . . 

If S is a field, then, by definition, x* = (1 — x)^ 1 for any i ^ 1. If S is 
an idempotent semiring, then, by definition 

x* = l®x®x 2 ®... = sup{l, x, x 2 , . . . , x n . . .}, (7) 

if this supremum (with respect to the canonical order =<;) exists. In this case, 
x* — 1 if x ^4 1. Therefore, x* = 1 in the semiring S^ 01 ' 6 ' for all x. For 
the semifield R max the closure operator x i— > x* is not defined for 1 -< x 
(however, R max can be supplemented by +oo, which turns this semifield into 
a semiring; in this case, x* = +oo for 1 -< x). It is clear that, for x =4 1, 
x* = 1 in R max , as well as in other idempotent semirings. These examples 
show that the closure x* of x is often calculated very simply for idempotent 
semirings. 

The closure operation for matrix semirings Mat n (S') can be defined and 
computed in terms of the closure operation for S; some methods are described 
in 0, ||, |6|, [IT|, |I2||. One such method is described below (LDM-factorization). 



The closure operation A A* in Mat n (S') satisfies identity (6), which implies 
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that if A* is defined, then X = A*B = A* © B is the solution to the matrix 
equation (1). 

Consider a nontrivial universal algorithm applicable to matrices over 
semirings with the closure operation defined. 

Example: Semiring LDM -Factorization 

Factorization of a matrix into the product A = LDM, where L and M are 
lower and upper triangular matrices with a unit diagonal, respectively, and D 
is a diagonal matrix, is used for solving matrix equations AX = B ||I3| . We 



construct a similar decomposition for the Bellman equation X = AX © B. 

For the case AX = B, the decomposition A = LDM induces the following 
decomposition of the initial equation: 

LZ = B, DY = Z, MX = Y. (8) 

Hence, we have 

A' 1 = M^D^L- 1 , (9) 

if A is invertible. In essence, it is sufficient to find the matrices L, D and M, 
since the linear system (8) is easily solved by a combination of the forward 
substitution for Z, the trivial inversion of a diagonal matrix for Y, and the 
back substitution for X. 

Using (8) CIS db pattern, we can write 

Z = LZ®B, Y = DY®Z, X = MX®Y. (10) 

Then 

A* = M*D*L*. (11) 

A triple (L, D, M) consisting of a lower triangular, diagonal, and upper tri- 
angular matrices is called an LD M- factorization of a matrix A if relations 
(10) and (11) are satisfied. We note that in this case, the principal diagonals 
of L and M are zero. 

The modification of the notion of LDM-factorization used in matrix anal- 
ysis for the equation AX = B is constructed by analogy with the construct 



suggested by Carre in [ 12 1 for L£7-factorization 
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We stress that the algorithm described below can be applied to matrix 
computations over any semiring under the condition that the unary opera- 
tion a i— > a* is applicable every time it is encountered in the computational 
process. Indeed, when constructing the algorithm, we use only the basic 
semiring operations of addition © and multiplication and the properties of 
associativity, commutativity of addition, and distributivity of multiplication 
over addition. 

If A is a symmetric matrix over a semiring with a commutative multipli- 
cation, the amount of computations can be halved, since M and L go into 
each other under transposition. 

We begin with the case of a triangular matrix A = L (or A = M) . Then, 
finding X is reduced to the forward (or back) substitution. 

Forward substitution 

We are given: 

• L = \\l l j\\i j =1 , where /* = for i < j (a lower triangular matrix with a 
zero diagonal); 

• b = m =1 . 

It is required to find the solution X = \\x' l \\^ =1 to the equation X = 
LX © B. The program fragment solving this problem is as follows. 

for i — 1 to n do 
{ x l := b l ; 

for j — 1 to % — 1 do 

x i ■= x i ® {1} @ x j ); } 

Back substitution 

We are given 

• M = ||?7i*-||"- =1 , where m* = for i > j (an upper triangular matrix 
with a zero diagonal); 
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It is required to find the solution X = ||a; l ||™ =1 to the equation X 
MX © B. The program fragment solving this problem is as follows. 

for % — n to 1 step —1 do 
{ x l := b' 1 ; 

for j = n to i + 1 step —1 do 

x % := x % © (m* 0i ! ); } 

Both algorithms require {n 2 — n)/2 operations © and 0. 

Closure of a diagonal matrix 

We are given 

• D = diag(di, . . . ,d n ); 

• b = pwu- 

It is required to find the solution X = \\x l \\f =1 to the equation X 
DX © B. The program fragment solving this problem is as follows. 

for % = 1 to n do 

x i := (di)* b l ; 

This algorithm requires n operations * and n multiplications 0. 

General case 

We are given 

• L — ||/j||™, =1 , where = if % < j; 

• D = diag(di, . . . ,d n ); 
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• M = \\nij\\2j = i, where rrij — if i > j; 

• B = pWtv 

It is required to find the solution X = \\x l \\f =1 to the equation X = 
AX(BB, where L, D, and M form the L.D M-factorization of A. The program 
fragment solving this problem is as follows. 

FORWARD SUBSTITUTION 

for % = 1 to n do 
{ x l := b l ; 

for j — 1 to i — 1 do 

x i ■= x* © (ij0jJ'); } 
CLOSURE OF A DIAGONAL MATRIX 
for i — 1 to n do 

x l := (di)* b l - 
BACK SUBSTITUTION 
for i = n to 1 step —1 do 
{ for j = n to i + 1 step —1 do 
:= x % © (m* x- 7 ); } 

Note that x % is not initialized in the course of the back substitution. The 
algorithm requires n 2 — n operations ©, n 2 operations 0, and n operations *. 

LDM-factorization 

We are given 

_ A — II n i lln 

It is required to find the LDM-factorization of A: L = ||^||"j =1 , D = 
diag(di, . . . ,d n ), and M = \\m l j\\2j =1 , where lj — if % < j, and = if 
i > j- 

The program uses the following internal variables: 

_ n — \\ n i \\n 

• ° - \\ c j\\i,j=i, 
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• d. 

INITIALISATION 

for i = 1 to n do 

for j = 1 to n do 
c) = a); 
MAIN LOOP 
for j = 1 to n do 
{ for i = 1 to j do 
f l := a*; 
for = 1 to j — 1 do 

for i = k + 1 to j do 

for i = 1 to j — 1 do 

a*- := (aj)* u*; 
aj := w- 5 ; 

for = 1 to j — 1 do 

for % = j ' + 1 to n do 

a} :=4©(4©t» fc ); 
d = (^')*; 

for i = j + 1 to n do 
a} := a* d; } 



This algorithm requires (2n 3 — 3n 2 + n)/6 operations ©, (2n 3 + 3n 2 — 
5n)/6 operations 0, and n(n + l)/2 operations *. After its completion, the 
matrices L, D, and M are contained, respectively, in the lower triangle, on 
the diagonal, and in the upper triangle of the matrix C. In the case when A 
is symmetric about the principal diagonal and the semiring over which the 
matrix is defined is commutative, the algorithm can be modified in such a 
way that the number of operations is reduced approximately by a factor of 
two. For details see [13| . 

[0- 



Other examples can be found in H], [11 
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6. SOFTWARE IMPLEMENTATION 
OF UNIVERSAL ALGORITHMS 



Object-oriented languages (e.g., C ++ and Java) and programming sys- 
tems that allow abstract data types to be defined provide convenient means 
for the software implementation of universal algorithms. In this case, pro- 
gram units can operate with abstract (and variable) operations and data 
types. Specific values of operations are determined by the input data types, 
these operations (and data types) are implemented by additional program 
units. Recently, this type of programming technique has been dubbed generic 
programming (see, e.g., [I], 0). To help automate the generic programming, 
the so-called Standard Template Library (STL) was developed in the frame- 
work of C ++ [El M. However, high-level tools, such as STL, possess both 
obvious advantages and some disadvantages and must be used with caution. 

Using the generic programming technique, a program package was devel- 
oped in C ++ for solving problems in linear algebra over fields and semirings 
(for various computer implementation of the corresponding numeric domains) 
and optimization problems on graphs. A hierarchy of abstract data types 
for basic numeric fields, rings, semifields, and semirings was developed for 
various computer representations. In particular, various versions of the ratio- 
nal arithmetic |7| can be used and computations can be performed with any 
given accuracy. Solving systems of linear Bellman equations over idempotent 
semirings (by various methods), standard optimization problems on graphs 
can be solved (the dynamic programming problem, shortest path problem, 
widest path problem, etc.), including interval versions of those problems [§], 
The system provides a basis for a more powerful program package based 
on universal algorithms |J . This system will be described in detail in subse- 
quent publications. 
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